///*
// * SystemState.cpp
// *
// *  Created on: 22 Aug 2011
// *      Author: Allan
// */
//
//#include "SystemState.h"
//
//// C++ includes
//#include <cmath>
//
//// GeoReact includes
//#include "ElectrolyteSolution.h"
//#include "Utilities/Algorithms.h"
//
//SystemState::SystemState()
//{}
//
//SystemState::SystemState(const Multiphase& multiphase, const ReactionSystem& eReactions) :
//System(multiphase, eReactions)
//{}
//	
//SystemState::SystemState(const Multiphase& multiphase, const ReactionSystem& eReactions, const ReactionSystem& kReactions) :
//System(multiphase, eReactions, kReactions)
//{}
//	
//SystemState::SystemState(const System& system) :
//System(system)
//{}
//
//SystemState::SystemState(const System& system, double T, double P, const VectorXd& n) :
//System(system)
//{
//	SetState(T, P, n);
//}
//
//void SystemState::SetState(double T, double P, const VectorXd& n)
//{
//	this->T = T;
//	this->P = P;
//	this->n = n;
//	
//	a = multiphase.Activities(T, P, n);
//	g = multiphase.ActivityCoefficients(T, P, n);
//}
//
//const double SystemState::Get(const string& quantity) const
//{
//	if(quantity == "IonicStrength" || quantity == "IonicStrength(Effective)")
//	{
//		const vector<string> aqueousSpecies = GetPhase("Aqueous").GetSpecies();
//		
//		ElectrolyteSolution sol(aqueousSpecies);
//		
//		return sol.EffectiveIonicStrength(GetPhaseComposition("Aqueous", n));
//	}
//	
//	if(quantity == "IonicStrength(Stoichiometric)")
//	{
//		const vector<string> aqueousSpecies = GetPhase("Aqueous").GetSpecies();
//		
//		ElectrolyteSolution sol(aqueousSpecies);
//		
//		return sol.StoichiometricIonicStrength(GetPhaseComposition("Aqueous", n));
//	}
//	
//	if(quantity == "pH")
//	{
//		double iH = multiphase["H[+]"];
//		
//		return -std::log10(a[iH]);
//	}
//	
//	cerr << "Error: The quantity type " << quantity << " is not valid." << endl; return 0.0;
//}
//
//const double SystemState::Get(const string& quantity, const string& species) const
//{
//	const Index ispecies = multiphase[species];
//	
//	if(quantity == "Moles" || quantity == "Concentration")
//		return n[ispecies];
//	
//	if(quantity == "Activity")
//		return a[ispecies];
//	
//	if(quantity == "ActivityCoefficient")
//		return g[ispecies];
//	
//	if(quantity == "PartialPressure")
//	{
//		const VectorXd ng = GetPhaseComposition("Gaseous", n);
//		
//		return P * (n[ispecies] / ng.sum());
//	}
//	
//	if(quantity == "SI" || quantity == "SaturationIndex")
//	{
//		cerr << "SystemState::Get(\"SI\") needs implementation." << endl; 
//
////		const Reaction& reaction = reactionMap.find(species)->second;
////		return std::log10(reaction.ReactionQuotient(a) / reaction.GetEquilibriumConstant()(T, P));
//	}
//	
//	if(quantity == "K" || quantity == "EquilibriumConstant")
//	{
//		const Reaction& reaction = reactionMap.find(species)->second;
//		
//		return reaction.GetEquilibriumConstant()(T, P);
//	}
//	
//	if(quantity == "KineticRate")
//	{
//		cerr << "SystemState::Get(\"KineticRate\") needs implementation." << endl;
//		
////		const Reaction& reaction = reactionMap.find(species)->second;
////		return reaction.KineticRate(T, P, n, a);
//	}
//	
//	cerr << "Error: The quantity type " << quantity << " is not valid." << endl; return 0.0;
//}
//
//void SystemState::Print() const
//{
//	System::Print();
//	
//	cout << "Temperature: " << T << " Celsius" << endl;
//	cout << "Pressure: " << P << " bar" << endl << endl;
//	
//	multiphase.PrintState(T, P, n);
//	
//	cout << endl << endl;
//	
//	eReactions.PrintState(T, P, n, a);
//	
//	cout << endl << endl;
//	
//	kReactions.PrintState(T, P, n, a);
//}
